Long Range Order at Low Temperature in Dipolar Spin Ice 
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Recently it has been suggested that long range magnetic dipolar interactions are responsible for spin 
ice behavior in the Ising pyrochlore magnets Dy2Ti207 and Ho2Ti207. We report here numerical 
results on the low temperature properties of the dipolar spin ice model, obtained via a new loop 
algorithm which greatly improves the dynamics at low temperature. We recover the previously 
reported missing entropy in this model, and find a first order transition to a long range ordered 
phase with zero total magnetization at very low temperature. We discuss the relevance of these 
results to Dy2Ti207 and Ho2Ti207. 
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In frustrated magnetism, the term spin ice was re- 
cently coined by Harris and coworkers 0| to describe the 
analogy that exists between the statistical physics of cer- 
tain geometrically frustrated Ising pyrochlore magnets, 
and proton ordering in the hexagonal phase of ice (Ih) 
|^|-Q. For the Ising pyrochlore systems Ho2Ti207 and 
Dy2Ti2 07, the Ho^+ and Dy^^^ rare earth magnetic mo- 
ments reside on a network of corner sharing tetrahedra 
(Fig. Ill) . Each moment is forced by single-ion anisotropy 
to lie along the axis joining the centers of the two tetra- 
hedra that it belongs to QH]. For a simple theoretical 
model considering only nearest neighbor ferromagnetic 
(FM) exchange, the groundstate is macroscopically de- 
generate, but is required to have two moments pointing 
in and two pointing out of every tetrahedron, a constraint 
that maps exactly the two short and long proton bonds 
and the ice-rules for their arrangement in Ih [^0. This 
nearest neighbor FM model shows no ordering and is 
characterized by a broad Schottky-like peak in the mag- 
netic specific heat m. 

Both Ho2Ti207 lla and Dy2Ti207 |§ show quahta- 
tive properties roughly consistent with the basic spin ice 
picture of the simple nearest neighbor FM model [^,0. 
However, it has been shown recently that rather than 
nearest neighbor FM exchange, it is surprisingly the large 
dipolar interaction present in these materials that is re- 
sponsible for their spin ice behavior ||q-p^. For a model 
which we call dipolar spin ice, with the long range nature 
of the dipolar interaction properly handled using Ewald 
summation techniques, numerical results show a lack of 
magnetic ordering down to very low temperatures p. 
Furthermore, the dipolar spin ice model agrees quanti- 
tatively very well with specific heat data for Dy2Ti2 07 
1^ and Ho2Ti207 [^, as well as neutron scattering mea- 
surements on the latter material |g]. In other words, 
while the simple nearest neighbor FM model provides a 
simple and qualitative understanding of the spin ice phe- 
nomenon, there is now strong evidence that the dipolar 



spin ice model with its long range dipolar interactions 
provides a quantitatively accurate description of experi- 
mental results on real materials 10,01. 




FIG. 1. The lower left 'downward' tetrahedron of the py- 
rochlore lattice shows Ising spins (arrows) whose axes meet 
at the middle of the tetrahedron. For clarity, black and white 
circles on the lattice points denote other spins. White repre- 
sents a spin pointing into a downward tetrahedron while black 
is the opposite. The entire lattice is shown in an ice-rules 
state (two black and two white sites for every tetrahedron). 
The hexagon (thick gray line) shows a minimal loop move, 
which corresponds to reversing all colors (spins) on the loop 
to produce a new ice-rules state. 

As in the case of Ih, for dipolar spin ice it is unclear 
whether the long range interaction should cause the ab- 
sence of ordering (and nonzero entropy) down to zero 
temperature. The dipolar interaction is itself FM at 
nearest neighbor, and is thus prone to spin ice correla- 
tions. However, a priori one might expect that its longer 
range component should lift the nearest neighbor degen- 
eracy and induce the selection of an ordered state within 
the ice-rules manifold. We show in this work that this 
is precisely the case. Specifically, the dipolar spin ice 
model with long range interactions does possess a unique 
groundstate (apart from trivial global symmetry oper- 



ations) which develops at very low temperature. How- 
ever, for local dynamical processes (such as single spin 
fluctuations), the development of this ground state is 
completely dynamically inhibited. As we discuss below, 
this occurs because of high energy barriers separating 
quasi-degenerate ice-rules states. This allows paramag- 
netic spin ice behavior to occur despite no special spin 
or space symmetry in the system which would a priori 
prevent magnetic ordering. In this paper we explore the 
low temperature ordering properties of dipolar spin ice 
by taking advantage of 'loop moves' incorporated into a 
standard Metropolis Monte Carlo algorithm, a method 
considered previously in the context of two-dimensional 
square ice models |13 . Such moves allow us to explore 



degeneracy lifting effects within the ice-rules manifold in 
an efficient manner, something which is not possible via 
single spin flip fluctuations. We present here strong nu- 
merical evidence for a first order phase transition at ex- 
tremely low temperature in the dipolar spin ice model in 
zero field that recovers the entire low temperature resid- 
ual magnetic entropy of the system. 

For the pyrochlore lattice with Ising spins defined by 
local axes, the Hamiltonian with nearest neighbor ex- 
change and long range dipolar interactions is IB-|lC 
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where the spin vector S^^ labels the Ising moment of mag- 
nitude 15*1 = 1 at lattice site i and local Ising [111] axis 
Zi discussed earlier. Here J represents the exchange en- 
ergy and D = (/xo/47r)(7^/x^/r^„. However, because of 
the local Ising axes, the effective nearest neighbor energy 
scales are Jnn = J/3 and D^n = 515/3. 

As described in Ref. M, the long range nature of the 
dipolar interactions can be handled conveniently by the 
Ewald method. In that work, extensive numerical anal- 
ysis via single spin flip Monte Carlo simulations found 
no evidence of a transition to long range order. Rather, 
short range order dominated by ice-rules correlations was 
observed down to low temperatures, similar to that found 
in the nearest neighbor FM model [ |l2[ . 

Qualitatively, the dynamics of both models appear to 
be very similar. As the temperature is lowered, signifi- 
cant thermal barriers are created by the energy cost in- 
volved in fluctuating out of the ice-rules manifold. With 
single spin flips, fluctuations between states within the 
ice-rules manifold are also reduced, as it is impossible to 
do so without first breaking the two-in/two-out ice-rules. 
Such thermal barriers produce non-trivial and extremely 
slow dynamics. If a unique groundstate exists within the 
plethora of ice-rules states (~ (3/2)^/^) of the dipolar 
spin ice model (Eq. 1), these thermal barriers make the 



probability of reaching it in a numerical simulation using 
conventional spin flips exceedingly small. Consequently, 
the question concerning the nature of the groundstate be- 
comes difficult to answer using standard numerical tech- 
niques, and a different procedure must be applied p3[ . 
Since we found in Ref. B that long range dipolar in- 
teractions give rise to spin ice behavior, we take as a 
starting point for identifying the low energy excitations 
(quasi zero modes) of Eq. (1) the exactly degenerate ice- 
rules state. This is entirely analogous to the approach 
taken in considering the so-called 'energetic ice models' 
in two-dimensional square ice models p3| . 
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FIG. 2. Specific heat data from simulations using single 
spin flips (squares) and combined single spin flips and loop 
moves (filled circles). Interaction parameters are Dy2Ti2 07 
values Jnn = -1.24 K and Dnn = 2.35 K (as in Ref. [9]). 
Inset: Finite size scaling of the specfic heat peak height as 
a function of system size L = 2, 3, 4, 5. The scaling behavior 
Cpoak(i) ~ a + bL^ is consistent with that expected for a first 
order phase transition. 

In Fig. ^we denote each site of the pyrochlore lattice by a 
white or black circle which represents a spin pointing into 
or out of a 'downward' facing tetrahedron, respectively. 
In this particular example, the spin configuration shown 
forms an ice-rules state that can be transformed into an- 
other ice-rules state by reversing all the colors (spins) on 
the loop denoted by the gray hexagon. In general, six 
spins form the shortest loop, while larger loops are also 
possible. A loop can be constructed by simply choosing 
a starting lattice site and tracing out a closed path that 
involves tetrahedra which have exactly two spins on the 
path (see Fig. 1). Furthermore, each pair of spins which 
are neighbors on the path are such that one is pointing 
into and the other pointing out of their shared tetrahe- 
dron. As seen in Fig. |l|, such a loop is constructed of 
alternating black and white circles. 

For our numerical study of the dipolar spin ice model, 
this type of 'loop move' was utilized in conjunction with 
conventional single spin flip dynamics. Speciflcally, such 
loops are identified by allowing a wandering path to form 



a loop whenever it encounters any previously visited site 
and ignoring any 'dangling' spins in the path. This allows 
for a large number of short loops to be created, with an 
average length that tends to a finite value as the system 
size is increased. As explained above for the dipolar sys- 
tem, such 'loop reversal' moves are not true zero modes, 
but involve a small gain or lowering of the energy (small 
compared to J„„ + Dnn) which is handled by a standard 
Metropolis algorithm |14|. 

Our numerical simulations for the dipolar spin ice 
model were carried out on system sizes up to 2000 spins 
(of cubic unit cell length L=5) with O(IO^) spin flips per 
spin and O(IO^) loop moves. For all interaction param- 
eters Jnn and Dnn which show spin ice behavior using 
single spin flip dynamics only (Jnn/-Dnn ^ —0.91) [g|, 
we find that the acceptance ratio of the loop moves in- 
creases at low temperature as the system enters the spin 
ice regime, before dropping to zero just below the tem- 
perature at which the system appears to undergo a very 
sharp first order phase transition to a long range ordered 
state obeying the ice-rules. 

In Fig. g we present specific heat data obtained for a 
system with interaction parameters Jnn and -Dnn identi- 
fied in Ref. Q for the spin ice material Dy2Ti207. Using 
a single spin flip Monte Carlo algorithm, spin ice correla- 
tions develop over a large temperature regime (signified 
by the broad peak around 1.1 K), before the system dy- 
namically slows down into a disordered ice-rules state 
at low temperature. Using the loop algorithm in com- 
bination with single spin flips, the higher temperature 
data is reproduced before a very sharp transition is ob- 
served at Tc ~ 0.18 K, with a latent heat observed at the 
transition. The energy probability distribution displays a 
double-peak feature in a narrow temperature region close 
to Tc, another indicator that the transition is first order. 
To assess in a more quantitative way the nature of the 
phase transition, a finite-size scaling study was done (see 
inset of Fig. 2). Because of the extremely sharp nature of 
the specific heat at Tc, the method of slowly cooling in a 
Monte Carlo simulation with discrete temperature steps 
could not give sufficiently accurate data to resolve Cpeak 
within reasonable computer time. To avoid this problem, 
simulations were performed in a multicanonical ensem- 
ble ||I^ at a single temperature near Tc. This data was 
then re-weighted using Ferrenberg and Swendsen's tech- 
nique [|l6| , which allowed us to obtain the appropriate 
thermodynamic quantities to any degree of temperature 
resolution required. 

The ordered phase is similar to that found in the or- 
der by disorder transition in the antifcrromagnctic FCC 
Ising model p7| . In that frustrated system, an ordering 
of antiferromagnetically stacked FM planes is found. For 
the dipolar spin ice system considered here, the ordering 
vector q lies parallel to one of the cubic axes directions, 
specifically q = (0, 0, 27r/a) or its starred directions. To 
construct the ordered state, first consider a starting tetra- 



hedron with its six possible ice-rules states. For a given 
ordering vector q, this tetrahedron selects one of the four 
possible spin configurations (two independent configura- 
tions and their spin-reversals. Si — > — S^), with a to- 
tal magnetic moment for the tetrahedron perpendicular 
to q. The entire ordered state may then be described 
by planes (perpendicular to q) of such tetrahedra. The 
wavelength defined by q physically corresponds to anti- 
ferromagnetically stacked planes of tetrahedra, where a 
given plane has tetrahedra of opposite configuration to 
the plane above and below it. In Fig. 3 we show one such 
groundstate with ordering vector q = (0, 0, 27r/a). 

The transition to such a groundstate structure can be 
characterized by the multi-component order parameter 
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Such a labeling is natural given that the pyrochlore lat- 
tice can be viewed as an FCC lattice with a 'downward' 
tetrahedral basis (see Fig. 1). Thus j labels the FCC 
lattice points of the pyrochlore lattice, and the index a 
sums over the four spins comprising the basis attached 
to each j. The index a labels the three possible sym- 
metry related q ordering vectors. For a given q^, as de- 
scribed above, there are two ice-rules configurations and 
their reversals which can each form a groundstate. Thus 
TO = 1, 2 labels these possibilities with the phase factors 
{0™} describing the given configuration m. Each Ising 
variable a^^i ^^as value 1 (-1) when a spin points into (out 
of) its downward tetrahedron labeled by j. 

As written in Eq.(2), ^t™ has six degenerate compo- 
nents, each of which can take on a value between and 1. 
Upon cooling through the transition, the system selects a 
unique ordered configuration, causing the corresponding 
component of '5™ to rise to unity and all others to fall to 
zero. The component selected by the ordering is equally 
likely to be any one of the six. Fig. 3 is a plot of (^} 

for two system sizes, where (^P) = Jj2^=i J2l=i (*")^ 
is the magnitude of the multi-component order parame- 
ter. For T < Tc the two lattice sizes produce identical 
order parameters. By contrast, (^) for the smaller lat- 
tice shows a somewhat more pronounced rounding near 
Tc, and an increased residual value for large T. These 
results show a clear discontinuity of the order parameter 
at Tc, and hence a first order transition to the long range 
ordered phase we have identified. 

For all values of Jnn /-Dnn within the dipolar spin ice 
regime P|, we find a low temperature phase transition 
to the state discussed above. The transition is driven 
by the long range dipolar interactions and, therefore, 
Tc ^ 0.1 8K is independent of the strength of the near- 
est neighbor exchange J [Tc/Dnn ^ 0.08). The obser- 
vation of a finite ordering temperature using the algo- 
rithm presented here demonstrates that long range dipo- 



lar interactions between Ising spins on the pyrochlore lat- 
tice have no special exact symmetry that allow a macro- 
scopically degenerate ground state. This conclusion is 
also suggested within a mean field analysis ^!^^ , which 
shows that as the truncation of long range dipolar inter- 
actions is pushed out to further distances (up to 10^ near- 
est neighbors), the maximal eigenvalues of the normal 
mode spectrum become only quasi- degenerate through- 
out the Brillouin zone, as opposed to the completely 
flat spectrum (and macroscopic degeneracy) we find for 
the nearest neighbor spin ice model [E0[. Furthermore, 
the quasi-degenerate eigenvalues of the mean field theory 
have a very weak dispersion which is maximal at the FCC 
zone boundary and, therefore, predicts the same ordering 
wavevector q found here. 
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FIG. 3. Temperature dependence of the order parameter 
{'^) defined above for system sizes L=3 (triangles) and L=4 
(squares). Inset: The q = (0, 0, 27r/a) groundstate projected 
down the z axis. The four tetrahedra making up the cubic unit 
cell 1 appear as dark grey squares. The light grey square does 
not represent a tetrahedron, however its diagonally opposing 
spins occur in the same plane. The component of each spin 
parallel to the z axis is indicated by a -f and — sign. 

The question remains as to what extent our conclu- 
sions apply to the real spin ice materials Ho2Ti207 ||] 
and Dy2Ti207 M. The dipolar spin ice model may be an 
accurate description of these materials even at extremely 
low temperatures, while it is also possible that perturba- 
tions, H' , exist beyond Eq. n^ which could induce another 
type of groundstate selection. Similar to fh however, ir- 
respective of the origin of any ordering, its actual obser- 
vation may depend critically on the dynamical behavior 
of the materials. The inability of single spin fluctuations 
to connect different ice-rules states in phase space shows 
that at low temperatures relaxation via local dynamics is 
extremely slow. For both Ho2Ti207 and Dy2Ti207, the 
transition temperature for the ordered phase observed in 
our simulations is well below the temperature at which 
single spin fluctuations over extended length scales (and 
out of the ice-rules manifold) are thermally frozen out. 
Thus, while theoretically an ordered phase induced by 
long range dipolar interactions between fsing spins on the 
pyrochlore lattice does exist, its experimental observation 
will depend acutely on the dynamical processes of the real 



materials. Furthermore, one requires that perturbations 
H' are negligible, ie. that H'/D^n < T^/Dnn < 0.08. 

In conclusion, we predict that in the dipolar spin ice 
model, which is in quantitative agreement with experi- 
mental data on real systems in the temperature regime 
investigated so far Mu, a very low temperature transi- 
tion to a zero total moment structure exists with recovery 
of all residual entropy. However, it is unlikely that such a 
phase can be arrived at via conventional local dynamics. 
These results suggest that Ising pyrochlore magnets with 
long range dipolar interactions provide an even deeper 
analogy with the proton ordering in hexagonal ice water 
Ih than previously suggested. 
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